################################################################################
class rap_states_lsj(rap_states):
  """ Store set of LSJ states. """
  def __init__(self,filename=None,**args):
    self.clear()
    if filename: self.loadfile(filename,**args)

# ------------------------------------------------------------------------------

  def _reprocess(self):
    """ Redefine senority-like tags. """
    if len(self) == 0: return
    indices=zip(*sorted(zip(self.v_energy.values(),self.v_energy.keys())))[1]
    allowedterms={}

    have={}
    tmp={}
    for idx in indices:
      (conf,term,sen,j2,par,en)=self.level(idx)

      # check for legal term
      if conf not in allowedterms: allowedterms[conf]=conf.allowed_terms()
      if term not in allowedterms[conf]:
        msg='invalid term for given configuration: %s %s'%(conf,term)
        raise ErrorRAP(msg)

      # check for legal level
      (g,l)=rlf_termdecode(term)
      if j2 < abs(2*l-g+1) or j2 > 2*l+g-1 or j2%2 != (2*l+g-1)%2:
        msg='invalid 2J value for term: %s %s 2J=%s'%(conf,term,j2)
        raise ErrorRAP(msg)

      # get and store seniority
      key=(conf,term)
      if key not in have: have[key]=0
      new=True
      for ij in range(1,have[key]+2):
        keyj=(conf,term,ij)
        if keyj not in tmp: break
        if j2 not in tmp[keyj]: 
          new=False
          break
      if new: 
        have[key]+=1
        tmp[keyj]=[]
      tmp[keyj].append(j2)
      self.v_seniority[idx]=ij
      if have[key] > allowedterms[conf][term]:
        msg='too many terms (%s) in configuration: %s %s; should have %s'%\
            (have[key],conf,term,allowedterms[conf][term])
        raise ErrorRAP(msg)

# ------------------------------------------------------------------------------

  def _parse_nist_line(self,line):
    """ 
    Function to extract level information (if any) from a NIST data file 
    (from webpage).  
    """
    if 'LIMIT' in line.upper():
      tmp=re.findall(r'(\d+\.\d+)',line)
      if not tmp: return False
      sym=line.split()[0].capitalize()
      zeff=None
      if sym in ATOMSYMBOLS: zeff=ATOMSYMBOLS.index(sym)
      ionth=float(tmp[-1])
      return (zeff,ionth)
    parts=line.split('|')
    if parts[0] == line: return False
    if parts[0].strip() == '':
      conf=None                  # None means no conf info
    else:
      tmp=re.findall(r'(\d+[a-z]\d*)',parts[0])
      if tmp == []: 
        conf=False               # False means can't understand conf
      else:
        conf=''
        for t in tmp: 
          l=re.findall(r'([a-z])',t)[0]
          l=SPECLETTERS.index(l.upper())
          tmp=re.findall(r'(\d+)',t)
          if len(tmp) == 1: t=t+'1'
          conf=conf+t+' '
        conf=conf[:-1]
    if parts[1].strip() == '':
      term=None                  # None means no term info
    else:
      tmp=re.findall(r'(\d+[A-Z])',parts[1])
      if tmp == []:
        term=False               # False means can't understand term
      else:
        term=tmp[0]
    tmp=re.findall(r'(\*)',parts[1])
    tmp=re.findall(r'(\d+)/2',parts[2])
    if tmp == []:
      tmp=re.findall(r'(\d+)',parts[2])
      if tmp == []: return False
      jay=float(tmp[0])
    else:
      jay=0.5*float(tmp[0])
    tmp=re.findall(r'(\d+\.\d+)',parts[3])
    energy=None
    if tmp != []: energy=float(tmp[0].replace('+x',''))
    j2=int(2.0*jay)
    return (conf,term,j2,energy)

# ------------------------------------------------------------------------------

  def clear(self):
    """ Erase all data. """
    rap_states.clear(self)
    self.v_cup='LSJ'
    self.v_conf={}
    self.v_term={}
    self.v_seniority={}
    self.v_j2={}

# ------------------------------------------------------------------------------

  def clear_levels(self):
    """ Just erase level data. """
    rap_states.clear_levels(self)
    self.v_cup='LSJ'
    self.v_conf={}
    self.v_term={}
    self.v_seniority={}
    self.v_j2={}

# ------------------------------------------------------------------------------

  def loadfile(self,filename,**args):
    """ Load in level information from level data file. """
    ft=rlf_filetype(filename)
    if ft == 'adf04':
      self.load_adf04(filename,**args)
    elif ft == 'olg':
      self.load_olg(filename,**args)
    elif ft == 'nist':
      self.load_nist(filename,**args)
    elif ft == 'chianti_elvlc':
      self.load_chianti(filename,**args)
    elif ft == 'x006':
      self.load_xstar(filename,**args)
    elif ft == 'atomdb_ev':
      self.load_atomdb(filename,**args)
    else:
      msg='do not know how to read file: %s'%filename
      raise ErrorRAP(msg)

# ------------------------------------------------------------------------------

  def reindex(self,remap=None):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  If remap is None, then reindex according
    to energy order.
    """
    imap=rap_states.reindex(self,remap)
    self.v_conf=dict([(imap[i],val) for (i,val) in self.v_conf.items()])
    self.v_term=dict([(imap[i],val) for (i,val) in self.v_term.items()])
    self.v_seniority=dict([(imap[i],val) for (i,val) in \
                          self.v_seniority.items()])
    self.v_j2=dict([(imap[i],val) for (i,val) in self.v_j2.items()])

# ------------------------------------------------------------------------------

  def configurations(self,levlist=None,imap=False):
    """ Return list of all configurations. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self._mapme,levlist)
    out=[]
    for i in levlist:
      if self.v_conf[i] in out: continue
      out.append(self.v_conf[i])
    return out

# ------------------------------------------------------------------------------

  def configuration(self,idx,imap=False):
    """ Return configuration of level with given index. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_conf[tidx]

# ------------------------------------------------------------------------------

  def terms(self,levlist=None):
    """ Return list of all terms. """
    if levlist is None: levlist=self.indices()
    out=[]
    for i in levlist:
      if self.v_term[i] in out: continue
      out.append(self.v_term[i])
    return out

# ------------------------------------------------------------------------------

  def term(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_term[tidx]

# ------------------------------------------------------------------------------

  def seniority(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_seniority[tidx]

# ------------------------------------------------------------------------------

  def j2range(self,levlist=None,imap=False):
    """ Return tuple containing the minimum and maximum 2J values. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self.mapme,levlist)
    j2s=[self.v_j2[i] for i in levlist if i is not None]
    return (min(j2s),max(j2s))

# ------------------------------------------------------------------------------

  def j2(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_j2[tidx]

# ------------------------------------------------------------------------------

  def jrange(self,levlist=None,imap=False):
    """ Return tuple containing the minimum and maximum 2J values. """
    (j2min,j2max)=self.j2range(levlist,imap)
    return (0.5*float(j2min),0.5*float(j2max))

# ------------------------------------------------------------------------------

  def j(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return 0.5*float(self.v_j2[tidx])

# ------------------------------------------------------------------------------

  def level(self,idx,imap=False):
    """ Return all information for single level (conf,term,alpha,j2,par,E) """
    tidx=self._mapme(idx,imap)
    if tidx is None: return (None,None,None,None,None,None)
    out=(self.v_conf[tidx],self.v_term[tidx],self.v_seniority[tidx],\
         self.v_j2[tidx],self.v_parity[tidx], self.v_energy[tidx])
    return out

# ------------------------------------------------------------------------------

  def label(self,idx,imap=False,spacer=None):
    """ Return level label; optional spacer can be used instead of space. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    j2=self.j2(idx,imap=imap)
    if j2%2 == 0:
      jstr='_'+str(j2/2)
    else:
      jstr='_'+str(j2)+'/2'
    out=str(self.v_conf[tidx])+' '+self.v_term[tidx]+jstr
    if self.v_seniority[tidx] > 1: out+='#'+str(self.v_seniority[tidx])
    if spacer is not None: out=out.replace(' ',spacer)
    return out

# ------------------------------------------------------------------------------

  def weight(self,idx,imap=False):
    """ Return weight of level. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.j2(idx,imap=imap)+1

# ------------------------------------------------------------------------------

  def add_level(self,conf,term,j2,energy,idx=False):
    """ Add new level to list """
    if idx is not False and idx in self:
      msg='index already exists: %s'%idx
      raise ErrorRAP(msg)
    tidx=idx
    if tidx is False: 
      tidx=1
      if len(self) > 0: tidx=max(self.indices())+1
    tconf=conf
    if not isinstance(tconf,rap_config): tconf=rap_config(tconf)
    if self.num_electrons() is None:
      self.v_nelec=tconf.num_electrons()
    else:
      if tconf.num_electrons() != self.num_electrons():
        msg='inconsistent number of electrons; expect/have: %s, %s'%\
            (self.num_electrons(),tconf.num_electrons())
        raise ErrorRAP(msg)
    (g,l)=rlf_termdecode(term)     # will raise error if bad term string
    if j2 < abs(2*l-g+1) or j2 > 2*l+g-1 or j2%2 != (2*l+g-1)%2:
      msg='invalid 2J value for term: %s 2J=%s'%(term,j2)
      raise ErrorRAP(msg)
    self.v_conf[tidx]=tconf
    self.v_term[tidx]=term
    self.v_seniority[tidx]=-1
    self.v_j2[tidx]=j2
    self.v_parity[tidx]=tconf.parity()
    self.v_energy[tidx]=energy
    self._reprocess()

# ------------------------------------------------------------------------------

  def add_levels(self,toadd,idxstart=None):
    """ 
    Add set of levels to structure. 

    The input list, toadd, contains tuples with: energy (Ry), index,
    configuration, term, and 2J.  If the indexing is to be automated, then set 
    index in the tuples to -1.
    """
    if len(toadd) == 0: return
    toadd.sort()
    idxcur=1
    if len(self) > 0: idxcur=max(self.indices())+1
    if idxstart is not None:
      if idxstart < idxcur: 
        msg='invalid starting index: %s'%idxstart
        raise ErrorRAP(msg)
      idxcur=idxstart
    tmp={}
    nechk=None
    for (en,idx,conf,term,j2) in toadd:
      tidx=idx
      if tidx in self or tidx in tmp: 
        msg='duplicate index: %s'%tidx
        raise ErrorRAP(msg)
      if tidx < 0:
        tidx=idxcur
        idxcur+=1
      tconf=conf
      if not isinstance(conf,rap_config): tconf=rap_config(conf)
      if nechk is None: nechk=tconf.num_electrons()
      if tconf.num_electrons() != nechk:
        msg='inconsistent number of electrons; expected/have: %s, %s'%\
            (nechk,tconf.num_electrons())
        raise ErrorRAP(msg)
      try:                       # raise error if bad term string
        (g,l)=rlf_termdecode(term)
      except ErrorRAP,msg:
        msg='problem with %s %s 2J=%s... %s'%(conf,term,j2,msg)
        warnings.warn(msg)
        continue
      if j2 not in range(abs(2*l-g+1),2*l+g,2):
        msg='invalid 2J value for term: %s %s 2J=%s'%(conf,term,j2)
        raise ErrorRAP(msg)
      tmp[tidx]=(en,tconf,term,j2)
    self.v_nelec=nechk
    for idx,(en,conf,term,j2) in tmp.iteritems():
      self.v_conf[idx]=conf
      self.v_term[idx]=term
      self.v_seniority[idx]=-1
      self.v_j2[idx]=j2
      self.v_parity[idx]=conf.parity()
      self.v_energy[idx]=en
    self._reprocess()

# ------------------------------------------------------------------------------

  def set_lmap(self,other,simple=False,minimize=False):
    """ Map to other rap_states_lsj instance. """
    if not isinstance(other,rap_states_lsj):
      msg='argument must be of type rap_states_lsj'
      raise ErrorRAP(msg)
    if self.num_electrons() != other.num_electrons():
      msg='inconsistent number of electrons; expected/have: %s, %s'%\
          (self.num_electrons(),other.num_electrons())
      raise ErrorRAP(msg)

    if simple:
      tmap={}
      for i in other:
        k=None
        for j in self:
          if self[j] == other[i]:
            k=j
            break
        tmap[i]=k
      self.v_lmap=tmap
      self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items()])
      return

    # group level indices by (conf,term,2J); seniority numbers not trusted
    fromlist={}
    for i in self:
      key=(self.configuration(i),self.term(i),self.j2(i))
      if key not in fromlist: fromlist[key]=[]
      fromlist[key].append(i)
    tolist={}
    for j in other:
      key=(other.configuration(j),other.term(j),other.j2(j))
      if key not in tolist: tolist[key]=[]
      tolist[key].append(j)

    if minimize:
      self.v_lmap=self._map_by_energy(other,fromlist,tolist)
    else:
      self.v_lmap=self._map_straight(other,fromlist,tolist)
    for j in other:
      if j not in self.v_lmap: self.v_lmap[j]=None
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items() if v is not None])
    for i in self:
      if i not in self.v_lmap_r: self.v_lmap_r[i]=None

# ------------------------------------------------------------------------------

  def cfilter(self,clist,subset=None,imap=False):
    """ Filter list of states based on configurations in clist. """
    if not isinstance(clist,list) and not isinstance(clist,tuple): clist=[clist]
    tclist=[]
    for c in clist:
      tc=c
      if not isinstance(tc,rap_config): tc=rap_config(tc)
      tclist.append(tc)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.configuration(i) in tclist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def tfilter(self,tlist,alpha=None,subset=None,imap=False):
    """ 
    Filter list of states based on terms in tlist; allow for optional 
    input for the seniority (alpha).
    """
    if not isinstance(tlist,list) and not isinstance(tlist,tuple): tlist=[tlist]
    if alpha is not None:
      if not isinstance(alpha,list) and not isinstance(alpha,tuple): 
        alpha=[alpha]

    ttlist=[]
    for t in tlist:
      tt=t
      if isinstance(tt,tuple): tt=rlf_termencode(tt[0],tt[1])
      ttlist.append(tt)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if alpha is not None and self.seniority(i) not in alpha: continue
      if self.term(i) in ttlist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def j2filter(self,j2list,subset=None,imap=False):
    """ Filter list of states based on 2J values in j2list. """
    if not isinstance(j2list,list) and not isinstance(j2list,tuple): 
      j2list=[j2list]
    tj2list=map(int,j2list)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.j2(i) in tj2list: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def jfilter(self,jlist,subset=None,imap=False):
    """ Filter list of states based on J values in j2list. """
    if not isinstance(jlist,list) and not isinstance(jlist,tuple): jlist=[jlist]
    tjlist=map(float,jlist)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.j(i) in tjlist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def symfilter(self,j2,p=None,subset=None,imap=False):
    """ Filter list of states based on symmetry. """
    if p is None:
      try:
        (jj2,pp)=j2
      except:
        msg='do not understand symmetry'
        raise ErrorRAP(msg)
    else:
      jj2=j2
      pp=p

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      tj2=self.j2(i)
      tp=self.parity(i)
      if tj2 != jj2 or tp != pp: continue
      out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def LS_average(self,omit_incomplete=True):
    """ 
    Average levels into LS terms contained in rap_states_ls object. 

    Indices will be set to energy-ordered.
    """
    agg={}
    for i in self:
      key=(self.v_conf[i],self.v_term[i],self.v_seniority[i])
      if key not in agg: agg[key]=(0.0,0.0)
      (esum,wgt)=agg[key]
      esum+=float(self.v_j2[i]+1)*self.v_energy[i]
      wgt+=float(self.v_j2[i]+1)
      agg[key]=(esum,wgt)

    emin=min([esum/wgt for (esum,wgt) in agg.values()])
    toadd=[]
    for (conf,term,alfa),(esum,wgt) in agg.iteritems():
      if omit_incomplete and rlf_termweight(term) > int(wgt): continue
      e=esum/wgt
      toadd.append((e-emin,-1,conf,term))

    out=rap_states_ls()
    out.set_nuclear_charge(self.nuclear_charge())
    out.set_core_energy(self.core_energy())
    out.set_ionization_threshold(self.ionization_threshold())
    out.add_levels(toadd)
    return out

# ------------------------------------------------------------------------------

  def search(self,conf=None,term=None,j2=None,parity=None,alpha=None,\
             energy=None):
    """ 
    Look for matching levels in structure based on configuration, term, and
    2J values.  If duplicate terms are present and alpha is not given, a list
    of all matching indices are returned.  If alpha is given, the index matching
    that seniority is returned.  If energy is given (and alpha is not), then
    the senority closest to that energy is returned.  If conf and term are 
    both None, then return matching index based on j2, parity, and energy order
    given by alpha.
    """
    if j2 is None: 
      if conf is None or term is None:
        msg="if no 2J is specified, must provide configuration and term"
        raise ErrorRAP(msg)
      if alpha is None: alpha=1
      levs=self.cfilter(conf)
      levs=self.tfilter(term,alpha,subset=levs)
      return levs

    levs=self.j2filter(j2)
    if conf is None or term is None:
      if conf is not None: parity=conf.parity()
      if parity is None: 
        raise RuntimeError('need parity when conf or term is None')
      levs=self.pfilter(parity,subset=levs)
      if levs == []: raise RuntimeError('matching level not found')
      if alpha is not None:
        return levs[alpha-1]
      elif energy is not None:
        tmp=[]
        for i in levs: tmp.append((abs(energy-self.energy(i)),i))
        return min(tmp)[1]
      return levs

    levs=self.cfilter(conf,subset=levs)
    levs=self.tfilter(term,subset=levs)
    if levs == []: raise RuntimeError('matching level not found')
    if len(levs) == 1: return levs[0]
    if alpha is not None:
      for i in levs: 
        if self.seniority(i) == alpha: return i
      raise RuntimeError('matching level not found')
    elif energy is not None:
      tmp=[]
      for i in levs: tmp.append((abs(energy-self.energy(i)),i))
      return min(tmp)[1]
    return levs

# ------------------------------------------------------------------------------

  def load_olg(self,filename):
    """ Load in level information from an AUTOSTRUCTURE olg file """
    fin=open(filename,'r')

    # first get closed subshells and number of configurations
    closed=[]
    while 1:
      line=fin.readline()
      if line == '': 
        msg='configuration list not found'
        raise ErrorRAP(msg)
      if 'CONFIGURATIONS' not in line: continue
      tmp=re.findall(r'IN(.+)CONFIGURATIONS',line)
      if tmp: break
    nconf=int(tmp[0])
    while 1:
      tmp=re.findall(r'(\d+)\s+(\d+)',line)
      if tmp == []: break
      closed.extend(tmp)
      line=fin.readline()

    # get list of all open subshells
    while 1:
      line=fin.readline()
      if line == '': 
        msg='could not find sub-shell list'
        raise ErrorRAP(msg)
      if 'INTERNAL ORBITALS' in line: break
    torbs=[]
    while 1:
      line=fin.readline()
      if 'CONFIGURATION' in line: break
      if 'N L' in line: 
        tmp=re.findall(r'(\d+)\s+(\d+)',line)
        if tmp: torbs.extend(tmp)
    orbs=[]
    for o in torbs:
      if o not in closed: orbs.append(o)
    base={}
    for o in orbs: base[o]=0

    ne=None
    conflist=[]
    for iconf in range(nconf):
      while 1:
        if line == '': 
          msg='configuration list not found'
          raise ErrorRAP(msg)
        if 'COMBINATIONS' in line:
          tmp=re.findall(r'COMBINATIONS(.+)',line)
          if tmp: break
        line=fin.readline()
      tmp=tmp[0]
      while 1:
        line=fin.readline()
        if 'SLATER' in line: break
        tmp+=line[:-1]
      tconf=base.copy()
      nl=re.findall(r'(\d+)\s+(\d+)',tmp)
      if not nl: 
        msg='could not understand %s configuration'%rlf_ordinal(iconf)
        raise ErrorRAP(msg)
      for ss in nl: tconf[ss]+=1
      aconf=''
      for o in orbs:
        ss=str(o[0])+SPECLETTERS[int(o[1])].lower()
        n=tconf[o]
        aconf+=ss+str(n)+' '
      conflist.append(rap_config(aconf,ne=ne))
      if ne is None: ne=conflist[0].num_electrons()

    # get nuclear charge
    while 1:
      if line == '':
        msg='file ended before retrieval of nuclear charge'
        raise ErrorRAP(msg)
      line=fin.readline()
      if 'ATOMIC NUMBER' in line: break
    tmp=re.findall(r'ATOMIC NUMBER\s*(\d+)',line)
    zeff=int(tmp[0])

    # get level information
    while 1:
      line=fin.readline()
      if line == '': 
        msg='no level data found'
        raise ErrorRAP(msg)
      if '(EK-E1)/RY' in line: break
    tmp=re.findall(r'(-?\d+\.\d+)\s*$',line)
    core_energy=None
    if tmp: core_energy=float(tmp[0])
    toadd=[]
    while 1:
      line=fin.readline()
      tmp=line.split()
      if tmp == []: break
      try:
        idx=int(tmp[0])
      except:
        break
      term=rlf_termencode(abs(int(tmp[4])),int(tmp[5]))
      j2=int(tmp[6])
      iconf=int(tmp[7])-1
      if iconf < 0: continue
      conf=conflist[iconf]
      energy=float(tmp[8])
      toadd.append((energy,idx,conf,term,j2))
    fin.close()
    self.set_nuclear_charge(zeff)
    self.set_core_energy(core_energy)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_adf04(self,filename):
    """ Load in level information from an ADAS adf04 file """
    fin=open(filename,'r')

    # basic info
    line=fin.readline()
    tmp=re.findall(r'([\d.]+)',line)
    z=int(tmp[1])
    q=int(tmp[2])
    ip=float(tmp[3])
    ne=z-q+1

    # read in level data
    isic=True
    isls=True
    isca=True
    toadd=[]
    while 1:
      line=fin.readline()
      if '-1' in line: break
      tmp=re.findall(r'(\d+)(.+)\((.+)\)(.+)\((.+)\)\s+(\d+[\.\d]*)',line)[0]
      idx=int(tmp[0])
      conf=rap_config(tmp[1],ne=ne)
      g=int(tmp[2])
      l=int(tmp[3])
      j2=int(2.*float(tmp[4]))
      e=float(tmp[5])/INVERSECM
      toadd.append((e,idx,conf,rlf_termencode(g,l),j2))

      cwgt=conf.weight()
      lwgt=g*(2*l+1)
      wgt=j2+1
      if wgt != cwgt: isca=False
      if wgt != lwgt: isls=False
      if wgt not in  range(abs(2*l-g+1)+1,2*l+g+1,2): isic=False
    fin.close()

    # check for correct coupling scheme
    if not isic:
      if isca: 
        msg='adf04 file seems to be configuration averaged'
        raise ErrorRAP(msg)
      elif isls:
        msg='adf04 file seems to be LS averaged'
        raise ErrorRAP(msg)
      else:
        msg='adf04 file not LSJ or a mistake in the J-values'
        raise ErrorRAP(msg)

    self.set_nuclear_charge(z)
    self.set_ionization_threshold(ip)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_nist(self,filename):
    """ Load in NIST data file (text file taken from website). """
    ne=None
    zeff=None
    ip=None
    toadd=[]
    ionth=[]
    fin=open(filename,'r')
    skipped=0
    for line in fin:
      tmp=self._parse_nist_line(line)
      if tmp:
        if len(tmp) == 2:
          if tmp[0] is not None: zeff=tmp[0]
          ionth.append(tmp[1])
        else:
          j2=tmp[2]
          en=tmp[3]
          if tmp[0] is False:
            conf=None
            term=None
            skipped+=1
            continue
          elif tmp[0] is None:
            if conf is None or tmp[1] is False: 
              conf=None
              term=None
              skipped+=1
              continue
            conf=toadd[-1][2]
            term=toadd[-1][3]
          else:
            if tmp[1] is False: 
              conf=None
              term=None
              skipped+=1
              continue
            conf=rap_config(tmp[0],ne=ne)
            term=tmp[1]
          if en == 0.0: ne=conf.num_electrons()
          toadd.append((en,-1,conf,term,j2))
    fin.close()

    # remove states with no energy
    for i in range(len(toadd)-1,-1,-1):
      if toadd[i][0] is None: toadd.pop(i)

    # any levels skipped
    if skipped > 0:
      plural=''
      if skipped > 1: plural='s'
      msg='%s level%s not understood and skipped'%(skipped,plural)
      warnings.warn(msg)

    if len(ionth) > 0: ip=min(ionth)
    self.set_nuclear_charge(zeff)
    self.set_ionization_threshold(ip)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_chianti(self,filename,priority="experimental",only=False):
    """ 
    Read energy levels from CHIANTI file.  CHIANTI energy level files contain
    both theoretical and experimental energies.  The priority option allows
    you to select which energy takes precedence (if available).  The only flag
    will restrict this function to storing only levels where the priority
    energy is present.
    """
    id_idx=0         # level index
    id_k=1           # configuration index
    id_conf=2        # configuration string
    id_g=3           # 2S+1
    id_lval=4        # L (integer)
    id_l=5           # L (string)
    id_j=6           # J
    id_wgt=7         # weight (2J+1)
    id_obcm=8        # observed energy (1/cm)
    id_obry=9        # observed energy (Ry); bad conversion, use 1/cm value
    id_thcm=10       # theoretical energy (1/cm)
    id_thry=11       # theoretical energy (Ry); bad conversion, use 1/cm value

    # attempt to get nuclear charge from filename
    tmp=re.findall(r'(\w+)_\d+\.elvlc',filename)
    if tmp:  
      zeff=ATOMSYMBOLS.index(tmp[0].capitalize())
      self.set_nuclear_charge(zeff)

    ne=None
    toadd=[]
    fin=open(filename,'r')
    for line in fin:
      tmp=line.split()
      if tmp[0] == '-1': break
      if len(tmp) != 12:
        raise RuntimeError('cannot undertand CHIANTI energy level file')
      idx=int(tmp[id_idx])
      conf=rap_config(tmp[id_conf],ne=ne)
      parity=conf.parity()
      g=int(tmp[id_g])
      l=int(tmp[id_lval])
      term=rlf_termencode(g,l)
      j2=int(tmp[id_wgt])-1
      ob=float(tmp[id_obcm])/INVERSECM
      th=float(tmp[id_thcm])/INVERSECM
      if priority[0] == 't':
        en=th
        if en == 0.0 and idx > 1: 
          if only is True: continue
          en=ob
      else:
        en=ob
        if en == 0.0 and idx > 1: 
          if only is True: continue
          en=th
      if ne is None: ne=conf.num_electrons()
      toadd.append((en,idx,conf,term,j2))
    fin.close()
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_xstar(self,filename):
    """ 
    Load level information from xstar type 006 file.
    """
    # column id's
    id_energy=7
    id_j2=8
    id_ip=10
    id_g=12
    id_l=13
    id_z=14
    id_ion=16
    id_id=15
    id_conf=17

    # read level list from file
    ion=None
    ne=None
    ldata=[]
    have=[]
    xid0=None
    fin=open(filename,'r')
    for line in fin:
      if line.strip() == '': continue
      data=line.split()

      confstr=data[id_conf-1]
      xid=int(data[id_id-1])
      ion=int(data[id_ion-1])
      g=abs(int(data[id_g-1]))
      l=int(data[id_l-1])
      j2=int(float(data[id_j2-1]))-1
      w=g*(2*l+1)
      if g == 1 or l == 0:
        cup='XX'
      elif w < j2+1:
        cup='CA'
        if 'CA' not in have: have.append('CA')
      elif w == j2+1:
        cup='LS'
        if 'LS' not in have: have.append('LS')
      else:
        j2min=abs(g-1-2*l)
        j2max=g-1+2*l
        okay=False
        for j2t in range(j2min,j2max+1,2):
          if j2t == j2: okay=True
        if okay:
          cup='IC'
          if 'IC' not in have: have.append('IC')
        else:
          cup='??'
      e=float(data[id_energy-1])
      if e == 0.0:
        ne=rap_config(confstr).num_electrons()
        zeff=int(data[id_z-1])
        ion_thresh=float(data[id_ip-1])
      chk=confstr.lower()
      if 'su' in chk or 'sp' in chk: continue 
      if 'co' in chk or 'cn' in chk: continue
      if 'fake' in chk: continue
      ldata.append([e/RYDBERG,cup,confstr,g,l,j2,xid])
    fin.close()
    if ne is None: raise RuntimeError('cannot find ground state')

    # id's for ldata
    id_e=0
    id_cup=1
    id_conf=2
    id_g=3
    id_l=4
    id_j2=5
    id_xid=6

    # get number of electrons from ground state and standardize configurations
    for i in range(len(ldata)):
      ldata[i][id_conf]=rap_config(ldata[i][id_conf],ne=ne)

    # get total weight of each configuration
    weight={}
    wtot={}
    cups={}
    for i in range(len(ldata)):
      conf=str(ldata[i][id_conf])
      if conf not in wtot: wtot[conf]=ldata[i][id_conf].weight()
      if conf not in cups: cups[conf]=[]
      if ldata[i][id_cup] not in cups[conf]: cups[conf].append(ldata[i][id_cup])
      if ldata[i][id_cup] in ['XX','LS']:
        if conf not in weight:
          weight[conf]=ldata[i][id_g]*(2*ldata[i][id_l]+1)
        else:
          weight[conf]+=ldata[i][id_g]*(2*ldata[i][id_l]+1)
      else:
        if conf not in weight:
          weight[conf]=(ldata[i][id_j2]+1)
        else:
          weight[conf]+=(ldata[i][id_j2]+1)

    # check if configurations have a common coupling scheme
    for c in cups.keys():
      tmp=cups[c][:]
      if 'XX' in tmp: tmp.remove('XX')
      if len(tmp) > 1:
        raise RuntimeError('mixed coupling scheme not supported')
      if tmp == [] and len(have) == 1: tmp=have[:]
      cups[c]=tmp[:]
      if cups[c] == []: cups[c]=['IC']

    # check that 'IC' configurations have all possible levels
    for c in weight.keys():
      if 'IC' not in cups[c]: continue
      elif weight[c] > wtot[c]:
        raise RuntimeError('configuration, '+str(c)+', has too many levels')

    # add levels having IC coupling
    toadd=[]
    for lev in ldata:
      idx=lev[id_xid]
      conf=lev[id_conf]
      sc=str(conf)
      if 'IC' not in cups[sc] and cups[sc] != ['XX']: continue
      term=rlf_termencode(lev[id_g],lev[id_l])
      j2=lev[id_j2]
      energy=lev[id_e]
      toadd.append((energy,idx,conf,term,j2))
    if len(toadd) == 0:
      msg='no LSJ data found in file, %s'%filename
      raise ErrorRAP(msg)
    self.set_nuclear_charge(zeff)
    self.set_ionization_threshold(ion_thresh/RYDBERG)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_atomdb(self,filename):
    """ 
    Load level information from ATOMDB FITS file.
    """
    # need pyfits
    rlf_have_pyfits(stop=True)

    a=pyfits.open(filename)
    adat=a[1].data
    header=a[1].header

    # check file type
    if header['HDUCLAS1'] != 'E_LEVEL':
      msg="invalid HDUCLAS1 in FITS file"
      raise ErrorRAP(msg)

    ionpot=None
    if 'ION_POT' in a[1].header:
      ionpot=float(a[1].header['ION_POT'])

    ionstr=header['ION_NAME']
    (z,ne)=rlf_iondecode(ionstr)
    q=z-ne

    # read in level data
    toadd=[]
    idx=0
    wgts={}                              # total configuration weight
    skipped=0
    for i in range(len(adat)):
      idx+=1
      g=int(2.*float(adat.field('S_QUAN')[i]))+1
      l=int(adat.field('L_QUAN')[i])
      j2=int(adat.field('LEV_DEG')[i])-1
      if g < 0:                                   # not LSJ
        skipped+=1
        continue
      e=float(adat.field('ENERGY')[i])            # in eV
      conf=rap_config(adat.field('ELEC_CONFIG')[i],ne=ne)

      if conf not in wgts: wgts[conf]=0
      wgts[conf]+=j2+1

      term=rlf_termencode(g,l)
      toadd.append((e/RYDBERG,idx,conf,term,j2))
    a.close()

    # check configuration weights
    for c in wgts:
      if c.weight() != wgts[c]:
        msg='incomplete configuration: %s'%str(c)
        warnings.warn(msg)

    if skipped > 0:
      msg='%d levels skipped; wrong coupling scheme'%skipped
      warnings.warn(msg)

    self.set_nuclear_charge(z)
    if ionpot is not None: self.set_ionization_threshold(ionpot)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def write_atomdb(self,filename,eref=None,comments=None):
    """ 
    Write level information to ATOMDB FITS file (from Adam Foster).
    """
    # need pyfits
    rlf_have_pyfits(stop=True)

    # prep comments
    if comments is None:
      comments=[]
    elif isinstance(comments,str):
      comments=[comments]
    else:
      comments=list(comments)

    # reference for energy levels
    if eref is None: eref=''

    # assemble energy level data
    elec_config=[]
    energy=[]
    e_error=[]
    n_quan=[]
    l_quan=[]
    s_quan=[]
    lev_deg=[]
    phot_type=[]
    phot_par=[]
    energy_ref=[]
    phot_ref=[]
    for ilev in self:
      (conf,term,alfa,j2,par,en)=self.level(ilev)
      valn=rlf_ssdecode(conf.valence_subshell())[0]
      (g,l)=rlf_termdecode(term)
      s=0.5*float(g-1)

      elec_config.append(str(conf))
      energy.append(RYDBERG*en)
      e_error.append(0.0)
      n_quan.append(valn)
      l_quan.append(l)
      s_quan.append(s)
      lev_deg.append(j2+1)
      phot_type.append(0)
      phot_par.append([0.0 for i in range(20)])
      energy_ref.append(eref)
      phot_ref.append('')

    # ion string
    z=self.nuclear_charge()
    ne=self.num_electrons()
    q=z-ne
    ionstr=rlf_ionencode(z,ne)

    # ionization potential
    ionpot=self.ionization_threshold()
    if ionpot is not None: ionpot*=RYDBERG

    #primary HDU, hdu0
    hdu0 = pyfits.PrimaryHDU()
    now = datetime.datetime.utcnow()

    hdu0.header.update('DATE', now.strftime('%d/%m/%y'))
    hdu0.header.update('E_LEVEL', "Atomic Energy Levels")
    hdu0.header.update('FILENAME', "Python routine")
    hdu0.header.update('ORIGIN', "ATOMDB",comment=os.environ['USER']+\
                       ", AtomDB project")
    hdu0.header.update('HDUCLASS', "ATOMIC",comment="Atomic Data")
    hdu0.header.update('HDUCLAS1', "E_LEVEL",comment="Atomic Energy levels")
    hdu0.header.update('HDUVERS', "1.1.0",comment="Version of datafile")

    #secondary HDU, hdu1:
    hdu1 = pyfits.new_table(pyfits.ColDefs(
    	      [pyfits.Column(name='ELEC_CONFIG',
    			     format='40A',
    			     array=elec_config),
    	       pyfits.Column(name='ENERGY',
    			     format='1E',
    			     unit='eV',
    			     array=energy),
    	       pyfits.Column(name='E_ERROR',
    			     format='1E',
    			     unit='eV',
    			     array=e_error),
    	       pyfits.Column(name='N_QUAN',
    			     format='1J',
    			     array=n_quan),
    	       pyfits.Column(name='L_QUAN',
    			     format='1J',
    			     array=l_quan),
    	       pyfits.Column(name='S_QUAN',
    			     format='1E',
    			     array=s_quan),
    	       pyfits.Column(name='LEV_DEG',
    			     format='1J',
    			     array=lev_deg),
    	       pyfits.Column(name='PHOT_TYPE',
    			     format='1J',
    			     array=phot_type),
    	       pyfits.Column(name='PHOT_PAR',
    			     format='20E',
    			     array=phot_par),
    	       pyfits.Column(name='ENERGY_REF',
    			     format='20A',
    			     array=energy_ref),
    	       pyfits.Column(name='PHOT_REF',
    			     format='20A',
    			     array=phot_ref)]
    	       ))

    hdu1.header.update('XTENSION', hdu1.header['XTENSION'],
    		    comment='Written by '+os.environ['USER']+now.strftime('%a %Y-%m-%d %H:%M:%S')+ 'UTC')
    hdu1.header.update('EXTNAME', ionstr, comment='Ion Name', before="TTYPE1")
    hdu1.header.update('HDUCLASS', 'ATOMIC',
    		    comment='Atomic Data', before="TTYPE1")
    hdu1.header.update('HDUCLAS1', 'E_LEVEL',
    		    comment='Energy level tables', before="TTYPE1")
    hdu1.header.update('ELEMENT', z,
    		    comment='Numer of protons in element', before="TTYPE1")
    hdu1.header.update('ION_STAT', q,
    		    comment='ion state (0 = neutral)', before="TTYPE1")
    hdu1.header.update('ION_NAME', ionstr, comment='Ion Name', before="TTYPE1")
    hdu1.header.update('N_LEVELS',len(energy) ,
    		     comment='Number of energy levels', before="TTYPE1")
    hdu1.header.update('HDUVERS1', '1.1.0',
    		     comment='Version of datafile', before="TTYPE1")
    if ionpot is not None:
      hdu1.header.update('ION_POT',ionpot,
             comment='Ionization potential (eV)',before="TTYPE1")

    for icmt in comments:
      hdu1.header.add_comment(icmt)

    # combine hdus
    hdulist = pyfits.HDUList([hdu0,hdu1])

    # write out file (overwrite any existing file)
    if os.path.isfile(filename): os.remove(filename)
    hdulist.writeto(filename, checksum=True)

# ------------------------------------------------------------------------------

  def write_xstar(self,filename,subset=None,precision=5):
    """ 
    Write xstar type 6 file. Must have ionization potential and nuclear charge
    specified.
    """
    if self.ionization_threshold() == None:
      raise ErrorRAP('must define ionization threshold')
    if self.nuclear_charge() == None:
      raise ErrorRAP('must define nuclear charge')
    indices=self.indices()
    if subset is not None: indices=subset

    zeff=self.nuclear_charge()
    ne=self.num_electrons()
    ionno=xstar_ionencode(zeff,ne)

    aform='%'+str(8+precision)+'.'+str(precision)+'E'
    iform='%6i'
    prefix=iform%6+iform%13+iform%0+iform%4+iform%6
    fout=open(filename,'w')
    for i in indices:
      conf=self.configuration(i)
      e=self.energy(i)*RYDBERG
      wgt=float(self.weight(i))
      ip=self.ionization_threshold()*RYDBERG
      np=-1.0
      if ip > 0.0: 
        z1=float(zeff-ne+1)
        de=abs(self.ionization_threshold()-self.energy(i))
        np=z1/sqrt(de)
      nval=0
      if conf is not None: nval=rlf_ssdecode(conf.valence_subshell())[0]
      (g,l)=rlf_termdecode(self.term(i))
      if g is None: 
        g,l=0,-1
      else:
        par=conf.parity()
        if par == 1: g=-g
      label=self.label(i,spacer='.')
      llen=len(label)
      line=prefix+iform%llen+' '+aform%e+aform%wgt+aform%np+aform%ip+' '+\
                  iform%nval+iform%g+iform%l+iform%zeff+iform%i+\
                  iform%ionno+' '+label+' %'
      fout.write(line+'\n')
    fout.close()

